home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / showpair.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  9.4 KB  |  494 lines  |  [TEXT/R*ch]
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include "clustalw.h"    
  6.  
  7.  
  8. /*
  9. *    Prototypes
  10. */
  11.  
  12. extern void *ckalloc(size_t);
  13. extern void  ckfree(void *);
  14. extern void error(char *,...);
  15. extern int *seqlen_array;
  16. extern char **seq_array;
  17.  
  18. void show_pair(void);
  19. static void make_p_ptrs(int *,int *,int,int);
  20. static void make_n_ptrs(int *,int *,int,int);
  21. static void put_frag(int,int,int,int);
  22. static int frag_rel_pos(int,int,int,int);
  23. static void pair_align(int,int,int);
  24. static void des_quick_sort(int *, int *, int);
  25.  
  26.  
  27. /*
  28. *     Global variables
  29. */
  30. extern int  dna_ktup, dna_window, dna_wind_gap, dna_signif; /* params for DNA */
  31. extern int prot_ktup,prot_window,prot_wind_gap,prot_signif; /* params for prots */
  32. extern int     nseqs;
  33. extern Boolean     dnaflag;
  34. extern double     **tmat;
  35. extern int     max_aa;
  36.  
  37. static int     next;
  38. static int     curr_frag,maxsf,vatend;
  39. static int     **accum;
  40. static int     *diag_index;
  41. static char     *slopes;
  42.  
  43. int ktup,window,wind_gap,signif;                  /* Pairwise aln. params */
  44. static int *displ;
  45. int *zza, *zzb, *zzc, *zzd;
  46.  
  47. extern Boolean percent;
  48.  
  49.  
  50. static void make_p_ptrs(int *tptr,int *pl,int naseq,int l)
  51. {
  52.     static int a[10];
  53.     int i,j,limit,code,flag;
  54.     char residue;
  55.     
  56.     for (i=1;i<=wind_gap;i++)
  57.            a[i] = (int) pow((double)(max_aa+1),(double)(i-1));
  58.  
  59.     limit = (int) pow((double)(max_aa+1),(double)ktup);
  60.     for(i=1;i<=limit;++i)
  61.         pl[i]=0;
  62.     for(i=1;i<=l;++i)
  63.         tptr[i]=0;
  64.     
  65.     for(i=1;i<=(l-ktup+1);++i) {
  66.         code=0;
  67.         flag=FALSE;
  68.         for(j=1;j<=ktup;++j) {
  69.             residue = seq_array[naseq][i+j-1];
  70.             if((residue<0) || (residue > max_aa)){
  71.                 flag=TRUE;
  72.                 break;
  73.             }
  74.             code += ((residue) * a[j]);
  75.         }
  76.         if(flag)
  77.             continue;
  78.         ++code;
  79.         if(pl[code]!=0)
  80.             tptr[i]=pl[code];
  81.         pl[code]=i;
  82.     }
  83. }
  84.  
  85.  
  86. static void make_n_ptrs(int *tptr,int *pl,int naseq,int len)
  87. {
  88.     static int pot[]={ 0, 1, 4, 16, 64, 256, 1024, 4096 };
  89.     int i,j,limit,code,flag;
  90.     char residue;
  91.     
  92.     limit = (int) pow((double)4,(double)ktup);
  93.     
  94.     for(i=1;i<=limit;++i)
  95.         pl[i]=0;
  96.     for(i=1;i<=len;++i)
  97.         tptr[i]=0;
  98.     
  99.     for(i=1;i<=len-ktup+1;++i) {
  100.         code=0;
  101.         flag=FALSE;
  102.         for(j=1;j<=ktup;++j) {
  103.             residue = seq_array[naseq][i+j-1];
  104.             if((residue<0) || (residue>max_aa)){
  105.                 flag=TRUE;
  106.                 break;
  107.             }
  108.             code += ((residue) * pot[j]);  /* DES */
  109.         }
  110.         if(flag)
  111.             continue;
  112.         ++code;
  113.         if(pl[code]!=0)
  114.             tptr[i]=pl[code];
  115.         pl[code]=i;
  116.     }
  117. }
  118.  
  119.  
  120. static void put_frag(int fs,int v1,int v2,int flen)
  121. {
  122.     int end;
  123.     
  124.     accum[0][curr_frag]=fs;
  125.     accum[1][curr_frag]=v1;
  126.     accum[2][curr_frag]=v2;
  127.     accum[3][curr_frag]=flen;
  128.     
  129.     if(!maxsf) {
  130.         maxsf=1;
  131.         accum[4][curr_frag]=0;
  132.         return;
  133.     }
  134.     
  135.         if(fs >= accum[0][maxsf]) {
  136.         accum[4][curr_frag]=maxsf;
  137.         maxsf=curr_frag;
  138.         return;
  139.     }
  140.     else {
  141.         next=maxsf;
  142.         while(TRUE) {
  143.             end=next;
  144.             next=accum[4][next];
  145.             if(fs>=accum[0][next])
  146.                 break;
  147.         }
  148.         accum[4][curr_frag]=next;
  149.         accum[4][end]=curr_frag;
  150.     }
  151. }
  152.  
  153.  
  154. static int frag_rel_pos(int a1,int b1,int a2,int b2)
  155. {
  156.     int ret;
  157.     
  158.     ret=FALSE;
  159.     if(a1-b1==a2-b2) {
  160.         if(a2<a1)
  161.             ret=TRUE;
  162.     }
  163.     else {
  164.         if(a2+ktup-1<a1 && b2+ktup-1<b1)
  165.             ret=TRUE;
  166.     }
  167.     return ret;
  168. }
  169.  
  170.  
  171. static void des_quick_sort(int *array1, int *array2, int array_size)
  172. /*  */
  173. /* Quicksort routine, adapted from chapter 4, page 115 of software tools */
  174. /* by Kernighan and Plauger, (1986) */
  175. /* Sort the elements of array1 and sort the */
  176. /* elements of array2 accordingly */
  177. /*  */
  178. {
  179.     int temp1, temp2;
  180.     int p, pivlin;
  181.     int i, j;
  182.     int lst[50], ust[50];       /* the maximum no. of elements must be*/
  183.                                 /* < log(base2) of 50 */
  184.  
  185.     lst[1] = 1;
  186.     ust[1] = array_size;
  187.     p = 1;
  188.  
  189.     while(p > 0) {
  190.         if(lst[p] >= ust[p])
  191.             p--;
  192.         else {
  193.             i = lst[p] - 1;
  194.             j = ust[p];
  195.             pivlin = array1[j];
  196.             while(i < j) {
  197.                 for(i=i+1; array1[i] < pivlin; i++)
  198.                     ;
  199.                 for(j=j-1; j > i; j--)
  200.                     if(array1[j] <= pivlin) break;
  201.                 if(i < j) {
  202.                     temp1     = array1[i];
  203.                     array1[i] = array1[j];
  204.                     array1[j] = temp1;
  205.                     
  206.                     temp2     = array2[i];
  207.                     array2[i] = array2[j];
  208.                     array2[j] = temp2;
  209.                 }
  210.             }
  211.             
  212.             j = ust[p];
  213.  
  214.             temp1     = array1[i];
  215.             array1[i] = array1[j];
  216.             array1[j] = temp1;
  217.  
  218.             temp2     = array2[i];
  219.             array2[i] = array2[j];
  220.             array2[j] = temp2;
  221.  
  222.             if(i-lst[p] < ust[p] - i) {
  223.                 lst[p+1] = lst[p];
  224.                 ust[p+1] = i - 1;
  225.                 lst[p]   = i + 1;
  226.             }
  227.             else {
  228.                 lst[p+1] = i + 1;
  229.                 ust[p+1] = ust[p];
  230.                 ust[p]   = i - 1;
  231.             }
  232.             p = p + 1;
  233.         }
  234.     }
  235.     return;
  236.  
  237. }
  238.  
  239.  
  240.  
  241.  
  242.  
  243. static void pair_align(int seq_no,int l1,int l2)
  244. {
  245.     int pot[8],i,j,k,l,m,flag,limit,pos,tl1,vn1,vn2,flen,osptr,fs;
  246.     int tv1,tv2,encrypt,subt1,subt2,rmndr;
  247.     char residue;
  248.     
  249.     if(dnaflag) {
  250.         for(i=1;i<=ktup;++i)
  251.             pot[i] = (int) pow((double)4,(double)(i-1));
  252.         limit = (int) pow((double)4,(double)ktup);
  253.     }
  254.     else {
  255.         for (i=1;i<=wind_gap;i++)
  256.                    pot[i] = (int) pow((double)(max_aa+1),(double)(i-1));
  257.         limit = (int) pow((double)(max_aa+1),(double)ktup);
  258.     }
  259.     
  260.     tl1 = (l1+l2)-1;
  261.     
  262.     for(i=1;i<=tl1;++i) {
  263.         slopes[i]=displ[i]=0;
  264.         diag_index[i] = i;
  265.     }
  266.     
  267.  
  268. /* increment diagonal score for each k_tuple match */
  269.  
  270.     for(i=1;i<=limit;++i) {
  271.         vn1=zzc[i];
  272.         while(TRUE) {
  273.             if(!vn1) break;
  274.             vn2=zzd[i];
  275.             while(vn2 != 0) {
  276.                 osptr=vn1-vn2+l2;
  277.                 ++displ[osptr];
  278.                 vn2=zzb[vn2];
  279.             }
  280.             vn1=zza[vn1];
  281.         }
  282.     }
  283.  
  284. /* choose the top SIGNIF diagonals */
  285.  
  286.     des_quick_sort(displ, diag_index, tl1);
  287.  
  288.     j = tl1 - signif + 1;
  289.     if(j < 1) j = 1;
  290.  
  291. /* flag all diagonals within WINDOW of a top diagonal */
  292.  
  293.     for(i=tl1; i>=j; i--) 
  294.         if(displ[i] > 0) {
  295.             pos = diag_index[i];
  296.             l = (1  >pos-window) ? 1   : pos-window;
  297.             m = (tl1<pos+window) ? tl1 : pos+window;
  298.             for(; l <= m; l++) 
  299.                 slopes[l] = 1;
  300.         }
  301.  
  302.     for(i=1; i<=tl1; i++)  displ[i] = 0;
  303.  
  304.     
  305.     curr_frag=maxsf=0;
  306.     
  307.     for(i=1;i<=(l1-ktup+1);++i) {
  308.         encrypt=flag=0;
  309.         for(j=1;j<=ktup;++j) {
  310.             residue = seq_array[seq_no][i+j-1];
  311.             if((residue<0) || (residue>max_aa)) {
  312.                 flag=TRUE;
  313.                 break;
  314.             }
  315.             encrypt += ((residue)*pot[j]);
  316.         }
  317.         if(flag) continue;
  318.         ++encrypt;
  319.     
  320.         vn2=zzd[encrypt];
  321.     
  322.         flag=FALSE;
  323.         while(TRUE) {
  324.             if(!vn2) {
  325.                 flag=TRUE;
  326.                 break;
  327.             }
  328.             osptr=i-vn2+l2;
  329.             if(slopes[osptr]!=1) {
  330.                 vn2=zzb[vn2];
  331.                 continue;
  332.             }
  333.             flen=0;
  334.             fs=ktup;
  335.             next=maxsf;
  336.         
  337.         
  338.         /*
  339.         * A-loop
  340.         */
  341.         
  342.             while(TRUE) {
  343.                 if(!next) {
  344.                     ++curr_frag;
  345.                     if(curr_frag>=FSIZE) {
  346.                         fprintf(stdout,"(Partial alignment)");
  347.                         vatend=1;
  348.                         return;
  349.                     }
  350.                     displ[osptr]=curr_frag;
  351.                     put_frag(fs,i,vn2,flen);
  352.                 }
  353.                 else {
  354.                     tv1=accum[1][next];
  355.                     tv2=accum[2][next];
  356.                     if(frag_rel_pos(i,vn2,tv1,tv2)) {
  357.                         if(i-vn2==accum[1][next]-accum[2][next]) {
  358.                             if(i>accum[1][next]+(ktup-1))
  359.                                 fs=accum[0][next]+ktup;
  360.                             else {
  361.                                 rmndr=i-accum[1][next];
  362.                                 fs=accum[0][next]+rmndr;
  363.                             }
  364.                             flen=next;
  365.                             next=0;
  366.                             continue;
  367.                         }
  368.                         else {
  369.                             if(displ[osptr]==0)
  370.                                 subt1=ktup;
  371.                             else {
  372.                                 if(i>accum[1][displ[osptr]]+(ktup-1))
  373.                                     subt1=accum[0][displ[osptr]]+ktup;
  374.                                 else {
  375.                                     rmndr=i-accum[1][displ[osptr]];
  376.                                     subt1=accum[0][displ[osptr]]+rmndr;
  377.                                 }
  378.                             }
  379.                             subt2=accum[0][next]-wind_gap+ktup;
  380.                             if(subt2>subt1) {
  381.                                 flen=next;
  382.                                 fs=subt2;
  383.                             }
  384.                             else {
  385.                                 flen=displ[osptr];
  386.                                 fs=subt1;
  387.                             }
  388.                             next=0;
  389.                             continue;
  390.                         }
  391.                     }
  392.                     else {
  393.                         next=accum[4][next];
  394.                         continue;
  395.                     }
  396.                 }
  397.                 break;
  398.             }
  399.         /*
  400.         * End of Aloop
  401.         */
  402.         
  403.             vn2=zzb[vn2];
  404.         }
  405.     }
  406.     vatend=0;
  407. }         
  408.  
  409. void show_pair(void)
  410. {
  411.     int i,j,k,dsr;
  412.     double calc_score;
  413.     
  414.     accum = (int **)ckalloc( 5*sizeof (int *) );
  415.     for (i=0;i<5;i++)
  416.         accum[i] = (int *) ckalloc(FSIZE * sizeof (int) );
  417.  
  418.     displ      = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
  419.     slopes     = (char *)ckalloc( (2*MAXLEN +1) * sizeof (char));
  420.     diag_index = (int *) ckalloc( (2*MAXLEN +1) * sizeof (int) );
  421.  
  422.     zza = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  423.     zzb = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  424.  
  425.     zzc = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  426.     zzd = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  427.  
  428.         if(dnaflag) {
  429.                 ktup     = dna_ktup;
  430.                 window   = dna_window;
  431.                 signif   = dna_signif;
  432.                 wind_gap = dna_wind_gap;
  433.         }
  434.         else {
  435.                 ktup     = prot_ktup;
  436.                 window   = prot_window;
  437.                 signif   = prot_signif;
  438.                 wind_gap = prot_wind_gap;
  439.         }
  440.  
  441.     fprintf(stdout,"\n\n");
  442.     
  443.     for(i=1;i<=nseqs;++i) {
  444.         if(dnaflag)
  445.             make_n_ptrs(zza,zzc,i,seqlen_array[i]);
  446.         else
  447.             make_p_ptrs(zza,zzc,i,seqlen_array[i]);
  448.         for(j=i+1;j<=nseqs;++j) {
  449.             if(dnaflag)
  450.                 make_n_ptrs(zzb,zzd,j,seqlen_array[j]);
  451.             else
  452.                 make_p_ptrs(zzb,zzd,j,seqlen_array[j]);
  453.             pair_align(i,seqlen_array[i],seqlen_array[j]);
  454.             if(!maxsf)
  455.                 calc_score=0.0;
  456.             else {
  457.                 calc_score=(double)accum[0][maxsf];
  458.                 if(percent) {
  459.                     dsr=(seqlen_array[i]<seqlen_array[j]) ?
  460.                             seqlen_array[i] : seqlen_array[j];
  461.                 calc_score = (calc_score/(double)dsr) * 100.0;
  462.                 }
  463.             }
  464. /*
  465.             tmat[i][j]=calc_score;
  466.             tmat[j][i]=calc_score;
  467. */
  468.  
  469.                         tmat[i][j] = (100.0 - calc_score)/100.0;
  470.                         tmat[j][i] = (100.0 - calc_score)/100.0;
  471.             if(calc_score>0.1) 
  472.                 fprintf(stdout,"Sequences (%d:%d) Aligned. Score: %lg\n",
  473.                            (pint)i,(pint)j,calc_score);
  474.             else
  475.                 fprintf(stdout,"Sequences (%d:%d) Not Aligned\n",
  476.                         (pint)i,(pint)j);
  477.         }
  478.     }
  479.  
  480.     for (i=0;i<5;i++)
  481.        ckfree((void *)accum[i]);
  482.     ckfree((void *)accum);
  483.  
  484.     ckfree((void *)displ);
  485.     ckfree((void *)slopes);
  486.     ckfree((void *)diag_index);
  487.  
  488.     ckfree((void *)zza);
  489.     ckfree((void *)zzb);
  490.     ckfree((void *)zzc);
  491.     ckfree((void *)zzd);
  492. }
  493.  
  494.